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ABSTRACT 

CN l/l/e study numerical methods of tomography in domains with a reflecting obstacle. It will be 

Q shown that tomography with sets containing both broken rays, i.e. rays reflecting at the ob- 

CN stacle, as well as unbroken rays, has a smaller error between the original and reconstructed 

^ image compared to classical tomography methods. 
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In biology, geophysics, oceanography, communications and other areas of science and tech- 



^ nology it is required to find the structure of a domain such as the human body, the Earth's 

2 interior, the oceans or the atmosphere from measurements taken at the domain's boundary. 

L_j This problem can be approached with the methods of tomography. When the domain contains 

P^ an obstacle new problems need to be addressed. In order to solve these problems, we study 

►> broken ray tomography. This work develops numerical methods using broken rays, or rays 

(^ reflecting at the obstacle, for solving the tomography problem. It will be shown that broken 

^ ray tomography leads to tomographic imaging with reduced error and better numerical stability 

1^ compared to classical tomography in the presence of reflecting obstacles. 
O 

^ Let u{x) designate acoustic wave propagation from point A to point B. Here u{x) denotes the 

LJ pressure at point x g M^. This process is governed by the wave equation 



utt — c(x)^An = 
where c{x) is the speed of sound at point x e M^. The travel time of wave propagation is 



nA,B)=f '' 



where the curve 71 is the acoustic geodesic from A to B. 

Assume that c(x) is close to a constant: c(x) = Co + e{x) where €{x) is small compared to 
Co- Then we will consider the approximation of 71 as the straight line segment 71. It is shown 
in (Romanov, 1974) that this linearization and more generally computing the travel time of 
u(x) over the known geodesic of the known speed Co can be justified under certain conditions. 



The proof given in (Romanov, 1974) is based on considering a perturbation parameter A and 
representing c{x) = Co + Xe{x). 

Let f{x) be a continuous function in Qq, wiiere Uq is a compact simply-connected set in M^ witii 
a smootii boundary. Consider dQo as tine observation boundary of Oq. Suppose we are given 
all integrals J f{l)dl = C, where 7 are all straight line segments or unbroken rays in Uq that 
have both of their endpoints in d^o. The classical Tomography Problem is to find f{x) in Qq 
knowing the values of all integrals C^. In the case of a domain Uq without an obstacle ili c ^0, 
this problem is widely studied theoretically and numerically (Natterer, 1986; Kak and Slaney, 
1988; Natterer and Wubbeling, 2001; Faridani, 2003). When there are obstacles present, 
this problem is much less studied. Some theoretical work is done in (Eskin, 2004; Eskin, 
2008; Eskin, 2006; Eskin, 2005) for domains with one obstacle and with both broken rays, i.e. 
rays reflecting at the obstacle, and unbroken rays. A key result from (Eskin, 2004) is that the 
tomography problem in the presence of one reflecting obstacle is well-posed. In this work, 
we implement a numerical method for solution of the tomography problem in a domain with 
one convex obstacle with piece-wise smooth boundary using unbroken and broken rays. If an 
obstacle Oi c ilo is present then we have the problem of recovering /(x) in r2o\ili with broken 
and unbroken rays 7. This problem is called the Broken Ray Tomography Problem(Eskin, 
2004). 

In a basic tomography setup transmitters and receivers of wave signals are placed at the do- 
main's boundary d^o- Signals are generated by the transmitters and received by the receivers. 
Travel times for signal propagation from transmitters to receivers are measured and, as dis- 
cussed, these travel times T{A,B) are the values of line integrals of a function f{x) = ^ 
where c{x) > is the speed of sound at point x g rio\f^i- This measurement procedure gives 
the C, data for solving the Tomography Problem by relating signal travel times to the values of 
line integrals of f{x). Given sufficient data, f{x) and from here the velocity c{x) are computed 
with tomographic reconstruction algorithms (Natterer, 1986; Kak and Slaney, 1988; Natterer 
and Wubbeling, 2001 ; Faridani, 2003). 

In nature and in technology, signals are often acoustic or electromagnetic waves. The propa- 
gation speed of electromagnetic waves is very high and this leads to very small travel times. In 
order to work with relatively larger travel times, I will start the description of the physical prob- 
lem by choosing high frequency acoustic waves or ultrasound as the signal carrier for obtaining 
tomography travel time data. The choice of ultrasound, or acoustic waves with frequency above 
20 kHz is crucial. The analysis and methods developed for ultrasound will be extended to other 
carriers and especially laser beams. The physical model for this work is linear ultrasonic wave 
propagation with reflection; a method for non-linear reconstructive ultrasound tomography is 
developed in (Schomberg, 1978). 

The acoustic geodesic for constant speed of sound is a straight line when there is no obstacle. 
Therefore, we consider two cases. In the first case, 7 = 71 is an unbroken ray. In the case of a 
broken ray, 7 = 71 U 72 is the union of two straight line segments that intersect at a reflection 



point at the obstacle. Reflection is mirror-like i.e. the angle of incidence is equal to the angle 
of reflection. 

For unbroken rays 

r^, . ^^ f ds f ds 

T{A,B)-- 



,^ Co + e{x{s)) J ^j^ Co + e{x{s)) 
and this model leads to the classical tomography problem with f{x) = -^-^^n^- 

For broken rays and a known obstacle, we know that the acoustic wave u{x) propagates along 
the known straight line segments 71 and 72. 

Then 

ds f ds f ds 



J^Co + e{x{s)) J^- Co + e{x{s)) L 



^ ^o + e{x{s)) J^^ Co + e{x{s)) J^^ Co + e(a:(s)) 

This model leads to the Broken Ray Tomography Problem. Reflection and the presence of 
an obstacle are unique to the broken ray model compared to the classical tomography prob- 
lem. The broken ray model is based on travel time computation and line integrals over broken 
and unbroken rays; the line integrals of classical tomography are over straight line segments. 
Moreover, the broken ray model requires tomographic reconstruction in the presence of an 
obstacle. For both models, we choose /(x) = ^^^^u\ ■ This physical model corresponds to a 
wave propagation environment that has almost constant wave propagation speed. Examples 
of such environments could be pollution plumes in the oceans or the atmosphere or regions of 
high concentration of minerals in water. Tomographic imaging of environmental pollution such 
as the recent oil spill in the Gulf of Mexico is a target application of broken ray tomography. 
The next chapter gives a geometric optics solution of the wave equation concentrated in a 
small neighborhood of a broken ray, which shows that we can work directly with broken rays 
when modeling mirror-like reflection of waves propagating approximately along broken rays. 

2 Geometric Optics Construction of a Brol^en Ray 

Following the approach in (Eskin, 2004; Eskin, 2008), I will construct geometric optics solutions 
of the wave equation, see (Eskin, 201 1), in a small neighborhood of a broken ray. 
Consider a wave, or signal, described by the wave equation 

wtt-c^AM = (2.1) 

where c = Co > is the constant speed of wave propagation and 

Ad<^x = 

for t > and Vt\{t) c ilo. where ili(t) is a moving convex obstacle with a piece-wise smooth 
boundary. We discussed in the introduction that under certain conditions such as small varia- 
tions of the speed of wave propagation, we can consider wave propagation along the geodesic 
of a wave with constant speed Co and that is why we consider a form of the wave equation with 
constant c. 



The geometric optics approacii to finding a solution of tiie wave equation tiiat approximates a 
broken ray is based on tine idea of looking for a solution of the form 

N N , 

p=0 p=0 

Starting from a plane wave solution e~*'='=*+*'="''^ of the wave equation, where k is large, i.e. 
A; -^ oo, and w = {wi,w2), \w\ = 1, is the unit direction vector of wave propagation, construct a 



geometric optics solution of 2.1 where 



N / . ^ 
-ikct+ikw-x V^ flpt^X,!, W) 

p=0 

approximates the broken ray before reflection as a sum of plane waves with varying amplitudes 
and 

kP 
p=0 



— ■tfi;cL-|-'i/i;■(/^^_■u;,x■y \ ^ 



approximates the broken ray after reflection as a sum of waves with a nonlinear phase tp{w, x). 
The term UN{x,t,w) makes the solution u{x,t) exact. 

Plugging the target solution in the wave equation leads to the equation 

p=0 p=0 

The amplitudes ap and bp are constructed for p = 0, 1,2, ... by plugging in the wave equation 
the target solution for A^ = 0, 1,2.... In the resulting expansion, equate to the sum of terms 
for equal powers of k for each of the two wave coefficients e"*^''*+*''"''^ and e-«*:ct+iA:^(w,x) 
This leads to a set of equations for the amplitudes ap{x,t,w) and bp{x,t,w). For example, for 
ao(x, t, w) and A^ = this step leads to the equation 

{^-c^A)ie-''^''+'''^--aoix,t,w)) = 

dt 

OXi^ 

+2ifct«ie-^^^*+^^'"-^^?^4?^^+ 

OXi 



dx2 

+2ife«;2e-*'='=*+'*^"'-"^?^4?^^) = 

0x2 



where w = {wi,w2), \w\ = 1, is the unit direction vector of wave propagation. Cancelling the 
term ^-''kct+tkw-x gj^gg 



~c\-kVa,ix,t,y.) + ^'°°J;;f'"^ + 2.fc^, ^«o(M,^) _ (2.4) 

2 2 ^ , N , 92ao(x,t,u;) 9ao(x,t,tf) 

-A; itJ2 ao(x,t, -u;) H — k h 2ikw2 ;^ ) = 

0x2 0x2 

Collecting terms and using the fact that w'^ = wi"^ + ti;2^ and that the coefficient for k in the 
above polynomial of k must be leads to 

dao{x,t,w) dao{x,t,w) _ 

dt dx 

Make the change of variables x = sw + tw±, where w± = {-W2,wi), and w -w^ =0. Then 

S = W • X 

T = 'W± • X 

Then the last equation can be written in the form 

dao{x,t,w) dao{x,t,w) _ 
dt ds 

The solution of this equation by the method of characteristics is 

ao{s, T, t, w) = f{s - ct, T, w) 
where / is an arbitrary smooth function. 

In order to localize the solution in a neighborhood of the broken ray, the solution is multiplied 
by cut-off functions. 



Continuing in this manner and plugging into the wave equation the sum 

^—ikct+ikw-x sr^N ap{x,t,w) -ikct+ikipiw ,x) s^N ^ \ TT ( ^ + »,A 

for successive values of A^ = 1, 2, ... and equating to in the resulting polynomial of k the sum 
of terms with a factor g-^kct+ikw-x fQ,. gqugi values of ^ leads to the recurrence relations 



d d d"^ 

2^c(^ +CW- —)aj{x,t,w) = (-"2 - c^A)aj-i{x,t,w) 



fori = 1,2,... 



Similarly for 6o(a;,t,tf) 

2^j -ikct+ikip(w,x) dbo[X, t, w) 

dt 

OX I 

, -ikct+ik'tl)(w,x) ^ bo[X,t,W) ^ 

dxi^ 

^.j^^^-ikct+iktbjw^x) dj^jx, w) dbojx, t, w) 
dxi dxi 

u2 ( d'4'\X., w) -^2^-ikct+ikip(w,x) u (^ + „„\ i ^-ikct+ikil){w,x) ^ bo{x,t,W) 

0x2 0x2^ 
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(2.5) 



r,„.j^^^-ikct+ik^(u,,x) dj'jx, w) dbojx, t, w) ^ ^ 
dx2 dx2 

Cancelling e~ikct+ik'^{w,x) ^ collecting terms and setting to the coefficients for the powers of k 
in the above polynomial of k leads to the eikonal equation 

IWP = 1 

and 

db()(x,t,w) dip(x,w) dbo(x,t,w) dil^(x,w) dbo(x,t,w) 

—^t '^^^, 9^^ + ^^^ a^^) = ° 

and in the general case for j > 1 the latter equation has the form 



;^M|^ - cV^ . Vb,ix,t,w)) = (|^ - c2A)6,_ 



In the recursive construction of a^ and bN, i.e. p = iV in the target solution sum, there are 
no terms to complete the recursive equations for the finite number of terms of order 0{^). 
Therefore, 

{^-c'A){UN{x,t,w)) = 0{^) 
and 

As discussed in (Eskin, 2005), this implies that Un = Oi-g^) 

The support of the constructed solution is in a small neighborhood of a broken ray. This neigh- 
borhood is defined by the cutoff functions in the solution. We can measure the travel time of 
an approximately linear signal sent from the beginning of the broken ray and received at the 
endpoint. Therefore, we can work directly with broken rays instead of their approximations. 



3 Numerical Methods for Broken Ray Tomography 

The broken ray tomography problem is well-posed(Eskin, 2004). This key result implies that 
numerical methods for image reconstruction through broken ray tomography should be more 
accurate and stable compared to the image reconstruction methods of classical tomogra- 
phy. Indeed, the classical tomography problem for a domain with an obstacle is ill-posed 
(Natterer, 1 986) and therefore, the well-posedness of the broken ray tomography problem is 
a crucial advantage for numerical methods for broken ray tomography. The well-posedness 
of the broken ray tomography problem follows from (Eskin, 2004) and the estimate from this 
paper 

/ |/(:,:)|^d. < C f° r(|2!^«£)lM>|. + |9^-(:.M.»W).,^,^ 

Jf^oVf^i Jo Jo OS d4> 

The notation for this estimate from (Eskin, 2004; Eskin, 2006) is as follows. Let 7^,0 be a 
broken or unbroken ray starting on point T in d^o with an endpoint R in r2o\r2i and direction 
angle 6 e S^ al point x = R. Let w{x,6) = J fds. Consider that the values of w{x,6) are 
known, and equal to C^, for all 9 e S^ when x e d^o- In this notation, 9{<f>) = (cos 0, sin 0), 
< < 27r, and /o is the length of d^o. 

This estimate implies that the solution of the broken ray tomography problem is uniformly 
bounded in the L^ norm. Indeed, the smooth function w{x,d) is defined for every x in the 
domain r2o\ili c Br, where Br = {x : \x\ < R}, of the function f{x) and < ^ < 27r. There- 
fore, its derivatives with respect to x and 9 are bounded. Therefore, ||/||L2(no\!^i) is bounded. 
Consider a perturbation of the input data wi{x,9) by a small amount to W2{x,9). Then the 
energy of the error can be estimated by 

i/,w - M.fd. < c I" r^fj^ummi _ *«Mf).Mi|.+ 

, , dwi{x{s), 9m dw2{x{s),9m ,2,^^^ 
+ '90 9,/. ' ' '*' 

The right hand side of the above inequality is bounded therefore the energy of the error 

/n„\ni 1/1 (^) - /2(a^)P is bounded. 

In developing numerical methods for solving the tomography problem in the presence of a 
reflecting obstacle we will use this result indirectly. Instead, we will directly rely on and build 
the numerical algorithms on the assumptions necessary for the above estimate and required 
by the theory of broken ray tomography. The key assumption and requirement of the theory of 
tomography in the presence of a reflecting obstacle is that all rays in the domain, both broken 
and unbroken, starting and ending at the observation boundary should be considered in order 
for the above estimate to hold. In other words, the problem is guaranteed to become well-posed 
when we consider all such rays. 



4 Broken Ray Tomography with the Kaczmarz Method 

In this work we apply the classical Kaczmarz methocl(Kaczmarz, 1937) to the Broken Ray 
Tomography Problem and show that broken ray tomography can be successfully performed 
with a well-known numerical method. We show through numerical experiments that linear 
systems corresponding to mixtures of broken and unbroken rays have more accurate solu- 
tions compared to linear systems corresponding to the same number of unbroken rays only. 
These numerical experiments show that broken ray tomography with the Kaczmarz method in 
the presence of a reflecting obstacle has much smaller error between the original and recon- 
structed image compared to classical tomography with the Kaczmarz method. 

It is well-known that reflections from an obstacle can improve the accuracy of tomographic 
imaging(Natterer, 201 0). The unique focus of the numerical solution of the Broken Ray Tomog- 
raphy problem is on determining the composition and properties of finite sets of rays that lead 
to linear systems that have more accurate and stable solutions when solved by the Kaczmarz 
method. We give a numerical solution of the broken ray tomography problem that restricts the 
solution of the tomography problem and equivalent problems to the solution of a class of linear 
systems that were obtained from sets of rays with favorable properties. 

I first describe the Kaczmarz method in the context of tomography in the presence of a reflect- 
ing obstacle and then I study the properties of finite sets of rays for solving the Broken Ray 
Tomography problem by the Kaczmarz method. The Kaczmarz method is the original and key 
method used in tomographic algebraic reconstruction algorithms and in tomography it is re- 
ferred to as ART. It is very flexible and works with rays with a wide range of known geometries. 

Consider a square domain M c M^ that contains the observable domain Qq i.e. $7i c f^o c 
M c M^. The square M is subdivided into a grid of N^ squares or cells of size d. We define 
f{x) = in Ui, inside the obstacle, and in M\i}o, outside the observation boundary, and look 
for a good approximation of f{x) in each cell of the grid. The value of f{x) is considered to be 
constant in each cell. We arrange linearly the N'^ cells into a column vector / = (/i^i, ..., /Ar,Af) 
where fij is the value of f{x) in cell M[i][j]. 

When a ray j intersects cell i of the vector f, then the length of the ray segment that the cell cuts 
from the ray is the weight of the cell with respect to this ray or t«[j][i]. The matrix of weights 
for all cells and all rays is denoted as W and has r rows and iV^ columns, where r is the total 
number of rays. Let Tj > be the travel time of ray j. Let 



j2^[mm = Tj 



This is the equation of a hyperplane in R^^ with normal vector Wj. Then the linear system of 
equations for all rays' travel-times can be expressed as 

Wf = T 



where T is the column vector of ray travel times. This linear system can be very large and even 
in initial numerical simulations has several thousand equations. The Kaczmarz method is an 
iterative method for solving large linear systems and it starts with an initial guess /(°). Then 

Wh -Wh 

where h = {i mod r) + 1 and r is the number of rays or number of rows of the linear system. 

Stefan Kaczmarz proved in (Kaczmarz, 1 937) the convergence of his iterative method for the 
solution of large regular linear systems. Tanabe proved in (Tanabe, 1971) that the Kaczmarz 
method will converge for any system of linear equations with nonzero rows. In general, esti- 
mating the Kaczmarz method's speed of convergence is still an open problem. 

The Kaczmarz method can be applied when the ray geometry is known and it is feasible to 
compute the intersection of the rays with the grid cells. In broken ray tomography the geometry 
of the broken ray is different from a straight line segment, however this geometry is known. 
In addition, the geometry and location of the obstacle are known and we can compute the 
intersections of the broken rays with the domain's cells outside the obstacle. Given an overde- 
termined system of rays with known geometries that includes both unbroken and broken rays, 
the Kaczmarz method can therefore be applied in order to reconstruct the velocity structure of 
an environment with a known obstacle. This is the proposed numerical solution of the broken 
ray tomography problem: 

Denote by A the set of all broken and unbroken rays that start and end at the observation 
boundary. Expand the set of rays that are used In the tomographic reconstruction to A 
or a discrete approximation of A. 

The theory of well-posedness of the broken ray tomography problem is based on consideration 
of the set of all rays, broken and unbroken, that start and end at the observation boundary. In 
other words, we need all rays if we want the reconstruction problem to be well-posed. There 
is an infinite number of such rays and for a numerical solution we need to approximate this 
condition in order to get as close as possible to the requirement for including all rays in the 
reconstruction. Our goal then is to characterize those finite sets of rays that approximate well 
the set of all rays. We start by defining what it means to approximate well. Consider an 
instance of the broken ray tomography problem with fixed domain Q.q, convex obstacle ili, 
function f defined in rio\f^i. and grid and cell size. Let s be a set of rays for reconstructing f 
that start and end at SOq. When signal travel time data is measured along the rays s, this set 
of rays leads to a linear system I. We solve I by the Kaczmarz method. Let / be the solution of 
I. Then 

e(^) = ll/-/ll 

is the error function of s induced by this instance of the broken ray tomography problem. In 
other words, e(s) is the error between the reconstructed and original value of f. The Euclidean 
and max norms or other norms can be chosen. 



Therefore, when solving a given instance of the broken ray tomography problem, we would like 
to select a set or sets s leading to a function e(s) that turns the given instance of the tomography 
problem into a well-posed one. The solution of the broken ray tomography problem found 
with the discrete set of rays s approximates the solution found with A. The theory of broken 
ray tomography implies that A leads to a well-posed instance of the broken ray tomography 
problem. 

Our goal is to have an efficient method of generating s or selecting from different si, S2, ... with- 
out solving the associated linear systems li,l2, ■■■■ Let two rays be equivalent if they intersect 
the same cells of the grid for a given instance of the broken ray tomography problem. This 
equivalence relation on the set of all broken and unbroken rays leads to a finite number of 
equivalence classes of rays. Indeed, each class is equivalent to a discrete ray composed from 
a finite number of cells from the grid. There are clearly a finite number of such discrete rays. 
This suggests a strategy for selecting s and a criteria for good approximation of A. Instead of 
seeking to minimize e(s), we select a finite set s such that each equivalence class of rays for 
the given instance of the broken ray tomography problem has at least one of its elements in s. 
Then s can be considered a finite set that approximates A. As a finite set of sets of measure 
0, s also has measure 0. For large grids, the number of elements in this set is a large but 
manageable instance of a discrete graph problem: 

Consider the cells of the grid that intersect d^o as the vertices of a graph d and the cells 
of the grid that intersect dQi as the vertices of a graph G2. A discrete unbroken ray is an 
edge between two vertices from d while a discrete broken ray is composed of two edges 
with one vertice from Gi and a shared vertice from G2. If the number of vertices of Gi is Vi 
and the number of vertices of G2 is V2 then the order of the number of elements of a discrete 
approximation of A is 

For a given instance of the broken ray tomography problem with an unknown f, it is possible to 
precompute such an approximation of A. This is possible because the geometry of all rays for 
a given instance of the broken ray tomography problem is known because of the linearization 
of ray propagation and the law of reflection. Data is collected for the rays in the set s and the 
resulting linear system solved by the Kaczmarz method. The justification for such a procedure 
is the improved accuracy and stability due to the well-posedness of the broken ray tomography 
problem. Such a finite set containing representatives from all equivalence classes of rays is 
still very large. In the numerical implementation described in the next chapter rays from this 
set are selected randomly without explicitly generating the whole set. 

5 Numerical Method Implementation 

Numerical methods for broken ray tomography in the presence of an obstacle should in the- 
ory be more accurate and stable reconstruction methods compared to classical tomographic 
reconstruction algorithms. This is an elegant result implied by the well-posedness of the bro- 
ken ray tomography problem proved in (Eskin, 2004). It implies that given a sufficiently large 



overdetermined linear system that corresponds to a large number of unbroken rays, we can im- 
prove the accuracy of the solution of the system by considering a large overdetermined system 
corresponding to the same number of rays a portion of which are broken. Moreover, broken 
rays are introduced naturally into the system due to reflection. 

Consider a fixed domain, grid, obstacle, signal wave carrier, and function / to be reconstructed. 
Consider all sets L of broken rays and unbroken rays in the domain. Some sets may contain 
only unbroken rays, some sets may contain only broken rays and other sets may contain both 
broken and unbroken rays. Each ray element of such a set is determined by its geometry and 
travel time between start and endpoint. Each such set determines a linear system for recon- 
structing f by tomographic reconstruction methods such as the Kaczmarz method. Therefore, 
to each such set of rays corresponds a linear system, or a set of equations, for reconstructing 
f. We conjecture that improved accuracy and stability of the reconstructed solution will become 
visible for a sufficiently large number of broken and unbroken rays that are uniformly distributed 
in the set of equivalence classes of all rays. This approach approximates the requirements of 
the theory of broken ray tomography which in turn implies well-posedness of the tomography 
problem in the presence of a reflecting obstacle. This section describes the implementation of 
our numerical solution of the broken ray tomography problem. 

We consider a uniform spatial distribution in which the broken and unbroken rays are uniformly 
distributed in the domain. Our goal is to approximate the set of all rays therefore other dis- 
tributions will favor some elements of the set of all rays and exclude other elements of this 
set. 

In order to verify the effectiveness of broken ray tomography with an obstacle, I generated an 
environment with a known obstacle and velocity that varies continuously in the observation 
region outside the obstacle. The following types of experiments were performed in order to 
compare broken ray tomography with tomography with rays that are straight line segments: a 
class of experiments for a fixed instance of the broken ray tomography problem and different 
ray sets approximating the set of all rays A, a class of exeriments where the size of the obstacle 
was varied, a class of experiments where the ratio of broken and unbroken rays was varied. 

A good numerical recipe for broken ray tomography should result in a reconstructed / that is 
closer in some norm to the true /. The basic algorithm for verifying the effectiveness of broken 
ray tomography can be summarized as follows: 

1 . Select a square domain M. 

2. Partition M into iV^ square cells of size d. In other words, the square grid M will have N cells 
of size d per row. 

3. Specify an observation boundary dQ.Q as a circle with center in M and configurable radius. 

4. Specify a known obstacle in M and represent it with the list of its boundary points. 

5. Generate a configurable number of transmitters and receivers located on the observation 
boundary. 

6. Generate broken rays by selecting randomly transmitters and points on the obstacle's 
boundary that are visible from the transmitters. Each pair of transmitter and obstacle boundary 



point can have exactly one corresponding receiver tiiat receives a signal sent from the trans- 
mitter that is reflected at the obstacle boundary point according to the law of reflection. In 
this work, we consider specular reflection i.e. reflection in a single direction such that angle of 
incidence is equal to the angle of reflection. We have performed and will report experiments 
with Lambertian reflection as well. 

In this network, there are no other blocking transmitters and receivers on the broken ray seg- 
ments connecting the reflection point on the obstacle with a transmitter and receiver pair. The 
choice of a circular observation boundary ensures this. 

Each transmitter receiver pair is considered only once: this corresponds to an implementation 
model in which at a given time a transmitter transmits in only one direction. In order to consider 
all rays, this restriction should be removed. 

7. Generate unbroken rays by selecting random pairs of transmitters and receivers that are 
not connected by broken rays: this corresponds to an implementation model in which at a 
given time a transmitter transmits in only one direction. In order to consider all rays, this 
restriction should be removed. Each pair has exactly one transmitter and exactly one receiver 
and the transmitter and receiver are visible from each other. For example, transmitters and 
receivers that are connected by a line segment that intersects the obstacle are not considered 
as endpoints of unbroken rays. 

8. Compute the weight matrix W of broken and unbroken ray intersections with the grid M. 
Elements w[j][i] of the matrix that are outside the observation boundary and inside the obstacle 
are set to 0. 

9. Choose a continuous function f{x,y) e M^. Discretize / by considering its domain as the 
cells of the grid. Set the value of / in each cell of the grid to be constant. The result /[I], 
f[2],...,f[N'^] are the values of / in the grid's cells. 

10. Use the known values of / from 9. to find the traveltime of a signal traveling along the 
broken or unbroken rays from 6. and 7. The physical model implies that the value of / in a 
given cell is equal to the reciprocal of the speed of ultrasound in the given cell for positive f 
or when the cell is inside the obstacle or outside the domain Qq. Therefore, multiplying the 
length of the section of the ray in the cell to the value of / in the cell gives the traveltime of 
ultrasound in the cell. By knowing the geometry of a ray, we know which cells of the grid the ray 
will intersect and we compute and add the traveltimes for each intersected cell to find the sum 
equal to the total traveltime of an ultrasonic signal along a given ray. This numerical integration 
along the path of ray j gives the traveltime T, of the ray for the chosen /. When this procedure 
is done for all rays, we know the vector T of travel times for all rays. 

11. Steps 8. and 10. give the coefficients matrix and right hand side of a linear system 
Wx = T. This system is very large and is solved with the Kaczmarz or other methods to find 
the value of the vector x. The vector x has the same meaning as the vector /, however, x is 
reconstructed numerically while / is chosen a priori. 

1 2. Now we can compare / and x and see how well the reconstructed value matches the true 
value of f. 

The numerical method is implemented in an original and custom Java software tomography 



framework called Euler. The program from (Gokgoz, Tanil and Onur, 2007) is another Java 
tomography framework for algebraic reconstruction that I have studied. 

6 Experimental Results 

Table[T|summarizes experimental results for tomographic reconstruction of a function f{x, y) = 
^\/{x - '^^Y + {y - yoY in a square domain with 4096 cells, or 64 cells per row where (xo, yo) 
is the center of the domain. The size of each cell is 13 points or units. The domain contains a 
square obstacle with 30 cells per row, and the circular observation boundary enclosing the ob- 
stacle has radius 350 units. There are 512 transmitters and 512 receivers along the boundary 
at equal angles between neighbor transmitters and the center of the observation boundary and 
between neighbor receivers and the center of the observation boundary. Each row of the table 
corresponds to a numerical solution for the function /(x, y) obtained with the same number of 
rays. Row 1 of TableQJshows results for an experiment with 1 26050 unbroken rays. These rays 
are an approximation of all possible unbroken rays between the transmitters and receivers that 
do not intersect the obstacle. Row 2 of Table [ijshows tomographic reconstruction for the same 
domain, obstacle and observation boundary with broken and unbroken rays. The total number 
of broken and unbroken rays in the domain is much larger compared to the total number of 
unbroken rays. In order to compare the effectiveness of broken ray tomography for the same 
number of rays, we take a random sample that has the same size 126050 as the number of 
unbroken rays in the solution with unbroken rays only. The ratio of broken and unbroken rays 
is 1:1 i.e. there are 63025 broken and 63025 unbroken rays in the sample that are chosen 
approximately uniformly. 

The error column shows the average error per cell between the original and reconstructed value 
of f. The number of iterations column shows the number of iterations of the Kaczmarz method 
before it converges to a solution. The criteria for convergence is that for at least two consecutive 
steps the currently computed value of f must be within a radius of convergence from the value 
of f computed at the previous step. I use the stronger max norm in addition to the classical 
Euclidean distance. In conclusion, the experimental results indicate that in the presence of an 
obstacle broken ray tomography is a more accurate approach to tomographic reconstruction 
compared to classical tomography. For example, visual inspection of the reconstructed images 
in figure |5] and figure |8] shows that compared to the image reconstructed with unbroken rays 
only, the reconstructed broken ray tomography image is significantly more accurate. 

Table[2]shows the performance of our numerical solution of the broken ray tomography problem 
and compares it to the performance of ART for a fixed instance of the tomography problem and 
ten different ray sets for each method. Reconstruction error increases when the fraction of 
unbroken rays in a ray set is close to 1 as shown by the results in Table[3]for the same instance 
of the broken ray tomography problem as in Table [3) Table |4]compares the performance of the 
two methods when the size of the obstacle is varied. 

In conclusion, we have observed in numerical experiments that broken ray tomography is on 
average three times more accurate compared to tomography without reflection. The theory of 



broken ray tomography predicts that the accuracy and advantages of broken ray tomography 
are much larger. 



Tomography Type 



ART 
BRT with Kaczmarz method 



Number of Rays 



Error 



Iterations 



126050 1 .80484955E-4 37144 

126050 4.820056689E-5 22728 



Table 1 : Comparison between ART and Broken Ray Tomography with the Kaczmarz method. 
We compare tomography without reflection and tomography with reflection in the presence of 
a reflecting obstacle. 



Experiment 


ART Error 


ART Iterations 


BRT Error 


BRT Iterations 


1 


1.104141e-004 


71502 


1.980839e-005 


77928 


2 


1.153358e-004 


69675 


9.159289e-005 


52365 


3 


1 .582274e-004 


46328 


2.191676e-005 


75798 


4 


9.90641 4e-005 


93348 


1. 62951 5e-005 


89842 


5 


1 .579209e-004 


45012 


2.107154e-005 


77348 


6 


1.511 201 e-004 


39476 


2.113383e-005 


80295 


7 


1 .349709e-004 


54259 


2.1 3951 7e-005 


78698 


8 


1.416369e-004 


57246 


1.001323e-004 


40882 


9 


1.313073e-004 


63075 


1 .732388e-005 


83562 


10 


1. 38371 9e-004 


52971 


2.1 89941 e-005 


75628 


Average 


1 .338370e-004 


59289.2 


3.525693e-005 


73234.6 



Table 2: Error and number of iterations for BRT for a fixed number of 126050 rays with 50% 
broken and 50% unbroken rays. The average error is 3.525693e-005 and the average number 
of iterations of the Kaczmarz method for finding a solution is 73234.600000. The results for 
ART tomography with 1 26050 unbroken rays are shown in the left two columns of the table. The 
average error for ART is 1 .338370e-004 and the average number of iterations is 59289.200000. 



Fraction of Unbroken Rays 


BRT Error 


BRT Iterations 


0.500000 


9.209899e-005 


58851 


0.550000 


1 .692470e-005 


84958 


0.600000 


2.651 309e-005 


60223 


0.650000 


5.372332e-005 


46297 


0.700000 


2.2991 21 e-005 


57166 


0.750000 


3.372089e-005 


41231 


0.800000 


2.778578e-005 


44129 


0.850000 


3.44301 6e-005 


33874 


0.900000 


4.696702e-005 


43283 


0.950000 


1. 286371 e-004 


47213 


Average 


4.837922e-005 


51722.5 



Table 3: Performance of broken ray tomography for a fixed instance of the tomography problem 
with ray sets of 1 26050 rays with different fractions of broken and unbroken rays. When the 
fraction of unbroken rays is close to 1 the reconstruction error increases. The average error is 
4.837922e-005 and the average number of iterations 51722.5. 



Side Length 


ART Error 


ART Iterations 


BRT Error 


BRT Iterations 


130 


1. 722081 e-004 


57490 


4.913426e-005 


83495 


156 


2.532327e-004 


41801 


4.810847e-004 


36541 


182 


2.786448e-004 


37833 


3.308767e-005 


85148 


208 


2.374432e-004 


55282 


2.977520e-005 


84985 


234 


1 .696454e-004 


135341 


3.608453e-005 


79498 


260 


2.21 861 3e-004 


76752 


2.544549e-005 


84010 


286 


1.888231 e-004 


91506 


2.249251 e-005 


84914 


312 


1.821 901 e-004 


79507 


2.4921 IOe-005 


81309 


338 


2.05381 5e-004 


51063 


1.41 6251 e-004 


31923 


364 


1 .743495e-004 


48376 


2.1 29221 e-005 


81740 


Average 


2.083780e-004 


67495.1 


8.649428e-005 


73356.3 



Table 4: Error and number of iterations for BRT for a fixed number of 126050 rays with 50% 
broken and 50% unbroken rays. Tomographic reconstruction is performed in ten experiments 
with different side lengths of the square obstacle. The average error is 8.649428e-005 and the 
average number of iterations of the Kaczmarz method for finding a solution is 73356.3. The 
results for ART tomography with 126050 unbroken rays are shown in the left two columns of 
the table. The average error for ART is 2.083780e-004 and the average number of iterations is 
67495.1. 




Figure 1 : Tomography with broken rays and unbroken rays 




Figure 2: Tomography with unbroken rays only 




Figure 3: Tomography with 126050 unbroken rays only 




Figure 4: Test case for tomography with 126050 unbroken rays only. Original value of f is 
continuously varying grey around the black square obstacle. 
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Figure 5: Tomography with 126050 unbroken rays only. Reconstructed f is inside the circular 
observation boundary and outside the black square obstacle. 




Figure 6: Tomography with 63025 broken and 63025 unbroken rays. 




Figure 7: Test case for tomography with 63025 broken and 63025 unbroken rays. Original 
value of f is continuously varying grey around the black square obstacle. 




Figure 8: Tomography with 63025 broken and 63025 unbroken rays. Reconstructed f is inside 
the circular observation boundary and outside the black square obstacle. 



7 Tomography in Three-Dimensional Domains with Obstacles 

The classical tomography problem for a three-dimensional domain without obstacles can be 
solved by taking two-dimensional plane cuts and solving the two-dimensional tomography prob- 
lem in each of the planar cuts. This approach is not generally applicable for broken ray tomog- 
raphy because for a given plane cut a ray that is in the plane of the cut before the reflection 
point is not guaranteed to be reflected in the same plane. Professor Eskin suggested a method 
for solving the three-dimensional tomography problem in a domain with one convex obstacle 
by considering plane cuts of Qq that do not intersect the obstacle. The tomography problem 
in each of these planes is well-posed and can be solved by the Radon transform and the al- 
gorithms of classical tomography (Natterer, 1986; Kaczmarz, 1937; Faridani, 2003; Natterer 
and Wubbeling, 2001). It is still interesting and practical to consider a solution of the three- 
dimensional tomography problem in a domain with an obstacle by plane cuts that intersect the 
obstacle. This approach is easier for applications because it is easier to choose the planes 
that cover the domain if we allow planes that intersect the obstacle. 

We can consider solving the three-dimensional broken ray tomography problem when the ob- 
stacle is a cylinder or a parallelipiped by taking parallel plane cuts that are perpendicular to the 
cylinder's axis or to one of the sides of the parallelipiped. All reflected rays will remain in the 
same plane as their respective incident rays. In the resulting plane cuts, the two dimensional 
obstacle will be a circle or a rectangle. Rectangles and circles have piece-wise smooth bound- 
aries required for obstacles by the theorems of broken ray tomography. It is an interesting 
question then whether generalized cylinders are the only three-dimensional shapes that ob- 
stacles can have so that the three-dimensional broken ray tomography problem can be solved 
by two-dimensional plane cuts which reduce the problem to the two-dimensional broken ray 
tomography problem in the union of the cuts that contains the whole obstacle. 

Proposition 7.1 . Define dO,i g R^ by F{x, y) = where F{x, y) is a smootin function tfiat is 
independent of z and suchi thiat ^i is convex and hias a smoothi boundary. Thie Broken Ray 
Tomographiy Problem in thie compact domain Oq c M^ and for an obstacle Oi c Uq can be 
solved by plane cuts of Oq a subset of which contains the obstacle ili . The cuts in that subset 
reduce the problem to the two dimensional broken ray tomography problem. 

Proof Let r^o c U Hz where n^ are planes parallel to the xy plane. Reconstruct the restriction 
of / to n^ and therefore f{x, y, z) in n^ by two-dimensional broken ray tomography in n^; when 
Uz does not intersect ^i f is determined by the Radon transform and can be reconstructed via 
the algorithms of classical tomography. Reconstruction by reduction to the two-dimensional 
broken ray tomography problem can be done because the surface normal at any point {xi,yi,z) 
on the boundary d^i is colinear with 

dF dF 

VF{xi,yi) = { — {xi,yi),—{xi,yi),0) 

and does not have a z component and is always contained in n^. Therefore, all incident 
rays in Qq f] U^ will be reflected in n^ because U^ contains the incident ray and surface normal 
VF(xi, yi) at the reflection point for that ray. Unbroken rays in n^ will remain in Uz by definition 



therefore Qq f] U^ contains all broken and unbroken rays that start on d^o f| n^ and that have 
their endpoints in d^o. n^fi'^^i is smooth and nzfl^i convex because dQi is smooth and 
ill convex and similarly II^ fi Qq is compact. Therefore / can be reconstructed by solving the 
well-posed two dimensional broken ray tomography problem in n^ fi Qq. 
The reconstructed restriction of / is unique in the plane n^ and the function f{x,y,z) is well- 
defined in Qq because the sets n^ are disjoint and each point {x,y,z) is in exactly one such 
set. D 

There are other shapes that allow reconstruction with plane cuts through the obstacle. For 
example, obstacle surfaces that are defined implicitly by functions that are independent of x or 
y can be cut by planes that are orthogonal to the x and y axis respectively. It is an interesting 
problem what all the shapes are for an obstacle so that the three dimensional broken ray to- 
mography problem can be solved by reducing it to the two dimensional broken ray tomography 
problem in planes intersecting the obstacle and the union of which contains the obstacle. 
One application of Proposition [Tj] is three dimensional tomographic reconstruction of the 
speed of sound in a domain that contains a cylindrical pipe. Imagine slices that cut the do- 
main and are perpendicular to the pipe. The thickness of each of these slices is small and they 
can be considered two dimensional. Moreover, the speed of sound in each of the slices can be 
considered close to a constant and this allows reconstruction by the numerical method from the 
previous chapter. The reconstructed values are combined and give f in the three dimensional 
region around the pipe. 
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